      subroutine tcmfdp1

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'cmfdp1.h'
      include 'sanm.h'
      include 'sanmjin.h'
      include 'setls.h'
      include 'eigv.h'
      include 'nem.h'
      include 'scmfdp1.h'

      real(NBF) jnetl1,jnetr1,jnetl2,jnetr2

! flux/current correction factor
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          do iem=1,2
            dtl(iem,j,im)=2*beta(iem,jm1,im)*beta(iem,j,im)/(beta(iem,jm1,im)+beta(iem,j,im))
            dtr(iem,j,im)=2*beta(iem,j,im)*beta(iem,jp1,im)/(beta(iem,j,im)+beta(iem,jp1,im))
          enddo          
        enddo
        dtl(1,1,im)=bcm_p1(1)
        dtr(1,tnode,im)=bcm_p1(2)
        dtl(2,1,im)=dtl(1,1,im)
        dtr(2,tnode,im)=dtr(1,tnode,im)
      enddo
      
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          btp2=beta(2,jp1,im)
          ppm1=psphi(1,jm1,im)
          pp1=psphi(1,j,im)
          ppp1=psphi(1,jp1,im)
          ppm2=psphi(2,jm1,im)
          pp2=psphi(2,j,im)
          ppp2=psphi(2,jp1,im)

          dhl(1,j,im)=-(jnet(1,j,im)+dtl(1,j,im)*(pp1-ppm1))/(pp1+ppm1)   ! current correction factor
          dhr(1,j,im)=-(jnet(1,jp1,im)+dtr(1,j,im)*(ppp1-pp1))/(ppp1+pp1)
          dhl(2,j,im)=-(jnet(2,j,im)+dtl(2,j,im)*(pp2-ppm2))/(pp2+ppm2)   
          dhr(2,j,im)=-(jnet(2,jp1,im)+dtr(2,j,im)*(ppp2-pp2))/(ppp2+pp2)

          sffdm1=(bt1*pp1+btm1*ppm1)/(bt1+btm1)  ! surface flux for the FDM
          sffdm2=(bt2*pp2+btm2*ppm2)/(bt2+btm2) 
          flcf(1,j,im)=(sflux(1,j,im)-sffdm1)/(pp1+ppm1)  ! flux correction factor
          flcf(2,j,im)=(sflux(2,j,im)-sffdm2)/(pp2+ppm2)
!          btm2=beta(2,jm1,im)
!          bt2=beta(2,j,im)
!          sffdm2=(bt2*avgflx(2,j,im)+btm2*avgflx(2,jm1,im))/(bt2+btm2)
          if (j.eq.tnode) then
            sffdm1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
            sffdm2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
            flcf(1,jp1,im)=(sflux(1,jp1,im)-sffdm1)/(pp1+ppp1)
            flcf(2,jp1,im)=(sflux(2,jp1,im)-sffdm2)/(pp2+ppp2)
          endif
        enddo
        if (bcb.eq.1) then
          dhl(1,1,im)=(0.375*sflux(2,1,im)-jnet(1,1,im))/sflux(1,1,im)-dtl(1,1,im)
          dhl(2,1,im)=(0.125*sflux(1,1,im)-jnet(2,1,im))/sflux(2,1,im)-dtl(2,1,im)
        endif
        if (bcu.eq.1) then
          dhr(1,tnode,im)=(-0.375*sflux(2,tnode+1,im)-jnet(1,tnode+1,im))/sflux(1,tnode+1,im) &
                         +dtr(1,tnode,im) 
          dhr(2,tnode,im)=(-0.125*sflux(1,tnode+1,im)-jnet(2,tnode+1,im))/sflux(2,tnode+1,im) &
                         +dtr(2,tnode,im) 
        endif
      enddo

! CMFD P1 for the 0th moment      
      do im=1,ng
        do j=1,tnode
          czf(1,j,im)=-dtl(1,j,im)+dhl(1,j,im)
          czf(2,j,im)=-dtr(1,j,im)-dhr(1,j,im)
          digf(j,im)=(dtr(1,j,im)-dhr(1,j,im)+dtl(1,j,im)+dhl(1,j,im))+xsr(j,im)*hz(j)
        enddo
      enddo

      do im=1,ng
         j=1
         del=digf(j,im)
         delinvf(j,im)=1/del
         jb=j
         do j=2,tnode
           deliauf(jb,im)=delinvf(jb,im)*czf(2,jb,im)
           del=digf(j,im)-czf(1,j,im)*deliauf(jb,im)
           delinvf(j,im)=1/del
           jb=j
         enddo
      enddo

!#define test
#ifdef test
      do im=1,ng
        do j=1,tnode
          tflx(j,im)=avgflx(1,j,im)
        enddo
        tflx(0,im)=avgflx(1,0,im)
        tflx(j,im)=avgflx(1,j,im)
      enddo
#endif

! source term
      do im=1,ng
        do j=1,tnode
          ss=0
          do m=1,ng         
            ss=ss+avgflx(1,j,m)*xssm(j,im)%from(m)
          enddo
!          src(1,j,im)=reigv*xchi(j,im)*psif(j)+ss*hz(j)
          src(1,j,im)=reigv*xchi(j,im)*psif(j)+(ss+2*xsr(j,im)*avgflx(2,j,im))*hz(j)
        enddo
!      enddo
!      do im=1,ng
! forward substitution
        j=1   
        rhs=src(1,j,im)
        y(j)=delinvf(j,im)*rhs
        jb=j
        do j=2,tnode
          rhs=src(1,j,im)-czf(1,j,im)*y(jb)
          y(j)=delinvf(j,im)*rhs
          jb=j
        enddo            
! backward substitution
        j=tnode   
        psphi(1,j,im)=y(j)
        jt=j 
        do j=tnode-1,1,-1
          psphi(1,j,im)=y(j)-deliauf(j,im)*psphi(1,jt,im)
          jt=j
        enddo
        do j=1,tnode
          avgflx(1,j,im)=psphi(1,j,im)-2*avgflx(2,j,im)
        enddo                  
      enddo ! for im          

      call dumflx

! CMFD P1 for the 2nd moment
      do im=1,ng
        do j=1,tnode
          czf(1,j,im)=-dtl(2,j,im)+dhl(2,j,im)
          czf(2,j,im)=-dtr(2,j,im)-dhr(2,j,im)
          digf(j,im)=(dtr(2,j,im)-dhr(2,j,im)+dtl(2,j,im)+dhl(2,j,im)) &
                    +(4*xsr(j,im)+5*xst(j,im))*de3*hz(j)
        enddo
      enddo

      do im=1,ng
         j=1
         del=digf(j,im)
         delinvf(j,im)=1/del
         jb=j
         do j=2,tnode
           deliauf(jb,im)=delinvf(jb,im)*czf(2,jb,im)
           del=digf(j,im)-czf(1,j,im)*deliauf(jb,im)
           delinvf(j,im)=1/del
           jb=j
         enddo
      enddo

! source term
      do im=1,ng
        do j=1,tnode
          ss=0
          do m=1,ng         
            ss=ss+avgflx(1,j,m)*xssm(j,im)%from(m)
!            ss=ss+tflx(j,m)*xssm(j,im)%from(m)            
          enddo
          src(1,j,im)=reigv*xchi(j,im)*psif(j)+(ss+2*xsr(j,im)*avgflx(2,j,im))*hz(j)
!          src(1,j,im)=reigv*xchi(j,im)*psif(j)+ss*hz(j)
          src(2,j,im)=-2*de3*(src(1,j,im)-xsr(j,im)*avgflx(1,j,im)*hz(j))
!          src(2,j,im)=-2*de3*(src(1,j,im)-xsr(j,im)*tflx(j,im)*hz(j))
        enddo
!      enddo
!      do im=1,ng
! forward substitution
        j=1   
        rhs=src(2,j,im)
        y(j)=delinvf(j,im)*rhs
        jb=j
        do j=2,tnode
          rhs=src(2,j,im)-czf(1,j,im)*y(jb)
          y(j)=delinvf(j,im)*rhs
          jb=j
        enddo            
! backward substitution
        j=tnode   
        psphi(2,j,im)=y(j)
        jt=j 
        do j=tnode-1,1,-1
          psphi(2,j,im)=y(j)-deliauf(j,im)*psphi(2,jt,im)
          jt=j
        enddo
        do j=1,tnode
          avgflx(2,j,im)=psphi(2,j,im)
        enddo                  
      enddo ! for im          

      call dumflx

! Determine the incoming partial current J_in
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          btp2=beta(2,jp1,im)
          ppm1=psphi(1,jm1,im)
          pp1=psphi(1,j,im)
          ppp1=psphi(1,jp1,im)
          ppm2=psphi(2,jm1,im)
          pp2=psphi(2,j,im)
          ppp2=psphi(2,jp1,im)
!          sffdml=(btm1*avgflx(1,jm1,im)+bt1*avgflx(1,j,im))/(btm1+bt1)
!          sffdmr=(btp1*avgflx(1,jp1,im)+bt1*avgflx(1,j,im))/(btp1+bt1)
!          sfluxl(1)=sffdml+flcor(j,im)*(avgflx(1,jm1,im)+avgflx(1,j,im))
!          sfluxr(1)=sffdmr+flcor(jp1,im)*(avgflx(1,jp1,im)+avgflx(1,j,im))  
!          sfluxl(2)=sflux(2,j,im)
!          sfluxr(2)=sflux(2,jp1,im)
!          sfluxl(2)=0
!          sfluxr(2)=0
!          jin(1,j,im)=0.25*sfluxl(1)+0.5*jnet_l+0.3125*sfluxl(2)
!          jin(2,j,im)=0.25*sfluxr(1)-0.5*jnet_r+0.3125*sfluxr(2)
!          avgflx(2,j,im)=0.8*(4*jin(1,j,im)-2*jnet(1,j,im)-sflux(1,j,im))
#define NEMP3
#ifdef NEMP3
! NEM P3 (using a pseudo flux)
          jnetl1=-dtl(1,j,im)*(pp1-ppm1)-dhl(1,j,im)*(pp1+ppm1)
          jnetr1=-dtr(1,j,im)*(ppp1-pp1)-dhr(1,j,im)*(ppp1+pp1)
          jnetl2=-dtl(2,j,im)*(pp2-ppm2)-dhl(2,j,im)*(pp2+ppm2)
          jnetr2=-dtr(2,j,im)*(ppp2-pp2)-dhr(2,j,im)*(ppp2+pp2)
          sffdml1=(btm1*ppm1+bt1*pp1)/(btm1+bt1)
          sffdmr1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
          sffdml2=(btm2*ppm2+bt2*pp2)/(btm2+bt2)
          sffdmr2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
          sfluxl(1)=sffdml1+flcf(1,j,im)*(ppm1+pp1)
          sfluxr(1)=sffdmr1+flcf(1,jp1,im)*(ppp1+pp1) 
          sfluxl(2)=sffdml2+flcf(2,j,im)*(ppm2+pp2)
          sfluxr(2)=sffdmr2+flcf(2,jp1,im)*(ppp2+pp2) 
!
          ajil(1,j,im)=0.25*sfluxl(1)+0.5*jnetl1-0.1875*sfluxl(2)
          ajir(1,j,im)=0.25*sfluxr(1)-0.5*jnetr1-0.1875*sfluxr(2)
          ajil(2,j,im)=0.4375*sfluxl(2)+0.5*jnetl2-0.0625*sfluxl(1)
          ajir(2,j,im)=0.4375*sfluxr(2)-0.5*jnetr2-0.0625*sfluxr(1)       
#endif
        enddo
        if (bcb.eq.1) jin(1,1,im)=0; ajil(1,1,im)=0; ajil(2,1,im)=0
        if (bcu.eq.1) jin(2,tnode,im)=0; ajir(1,tnode,im)=0; ajir(2,tnode,im)=0  
      enddo
 
      end subroutine
